CUBATURE ON WIENER SPACE IN INFINITE DIMENSION 
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Abstract. We prove a stochastic Taylor expansion for SPDEs and apply this 
result to obtain cubature methods, i.e. high order weak approximation schemes 
for SPDEs, in the spirit of T. Lyons and N. Victoir. We can prove a high-order 
weak convergence for well-defined classes of test functions if the process starts 
at sufficiently regular initial values. We can also derive analogous results in 
the presence of Levy processes of finite type, here the results seem to be new 
even in finite dimension. Several numerical examples are added. 



1. Introduction 

Let H he a, real separable Hilbert space. We consider the stochastic (partial) 
differential equation for a diffusion process X with values in H 

d 

(1) dX^ = (AXf + a{X-^))dt + f^^mdBl 

1=1 

or ~ in the presence of jumps - the stochastic differential equation for a jump- 
diffusion process X with values in H 

d e 

(2) dX^ = {AX^^ + a{X^. ))dt + ^ )dBl -f ^ 5, )dLl , 

i=i j=i 

where, in general, A : C H H denotes an unbounded linear operator, 

a, Pi, . . . , f]d,Si, . . . ,Se ■ H ^ H denote C°°-bounded vector fields - i.e. the vec- 
tor fields are smooth and all the derivatives (of degree > I) are bounded - and 
{Bt)t>o — {B} , ■ ■ ■ ,Bf)t>o denotes a finite dimensional Brownian motion on the 
Wiener space (i7, J^, P), and {Ll)t>o a compound Poisson process with jump-rate 
fij for j — 1, . . . , e. The initial value x G H appears as superscript in the notation 
X^ of the solution process, i.e. Xq = x. We assume that A is the generator of a 
Co-semigroup denoted by {St)t>o- The main reference for equations of the type ([T|) 
is the monograph of G. da Prato and J. Zabczyk [S]. In the (general) Levy case 
we refer to [6] and [19] , even though we do not need results of their strength in our 
case since all Levy processes under consideration are of finite type. 

Due to the unboundedness of the operator A several concepts of solutions to ^ 
arise. We refer for the precise definitions to the excellent monograph [THj. The 
most direct analogue to the finite dimensional setting is the concept of a strong 
solution, which is defined by 



(3) Xr = x+ 



/ {AX:^+a{X:^))ds + J2 f m!-)dBl + 
Jo Jo 

+ E fs,{X^-)dLl, 

.7 = 1 
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i.e. by the integrated version of equation Besides the obvious integrability 

assumptions we need in particular that e for alH > almost surely with 

respect to P. More relevant is the concept of a mild solution, which is related to 
^ via variation of constants. Indeed, a mild solution is a process satisfying 



given the obvious integrability assumptions. By Ito's formula, any strong solution 
is a mild solution, but not the other way round. Note that mild solutions need 
not be semi-martingales, because there is no semi-martingale decomposition if the 
process evolves outside of 'D{A). Consequently, a Stratonovich formulation of ^ 
is, in general, not possible for mild solutions. 

Of course, neither strong nor mild solutions can usually be given explicitly, which 
makes numerical approximation necessary. We are interested in weak approximation 
of the solution in the sense that we want to approximate the value 



for a suitable class of test functions f : H ^ M. at initial values x ^ H. It is 
well-known that the function {t,x) ^ Ptf{x) solves the Kolmogorov equation in 
the weak sense, see for instance [5] in the diffusion case. 

Let us assume for a moment that there are no jump-components: usually, infi- 
nite dimensional SDEs are numerically solved by finite element or finite difference 
schemes, see, for instance, [H], [23], [S] and [TUj. For HJM models of financial 
mathematics a finite difference scheme and a finite element scheme have been im- 
plemented in [3]. This means that the original equation is projected onto some 
finite dimensional subspace Hh C H and A is approximated by some operator Ah 
defined thereon. This procedure, which corresponds to a space discretization of the 
stochastic PDE, is followed by a conversion of the stochastic differential equation 
on Hh to a stochastic difference equation by discretizing in time, using an Euler 
method or a related scheme. Finally, the stochastic difference equation is solved 
by Monte-Carlo simulation, which may be interpreted as a discretization on the 
Wiener space. For general information about approximation of finite dimensional 
SDEs see [H]. 

We want to tackle the problem in the reverse order: we want to do the discretiza- 
tion on the Wiener space first, reducing the problem to a deterministic problem, 
i.e., one replaces the d-dimensional Brownian motion with finitely many trajecto- 
ries of bounded variation chosen with well-defined probabilities such that certain 
moments (of iterated Stratonovich integrals) match. The resulting deterministic 
problem can be solved by standard methods for the numerical treatment of deter- 
ministic PDEs, e.g. by standard finite element or finite difference methods. The 
benefit is that once the discretization on the Wiener space has been done, we can 
use the well-established theory of the corresponding deterministic PDE-problems, 
without any complications from stochasticity. Our method of choice for discretiza- 
tion on is ^'cubature on Wiener space", developed by Terry Lyons and Nicolas 
Victoir in [TH] and by Shigeo Kusuoka in [H] and [H] , see also [12] . In the spirit of 
these methods we shall obtain weak approximation schemes of any prescribed order 
of convergence. Notice here that we discretize in the presence of the unbounded 
operator A in the drift vector field. Certainly our Assumptions 12.11 seem very re- 
strictive, but these assumptions are the appropriate analogues of the assumptions 
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in finite dimension that the vector fields are bounded, C°°-bounded (see [TJ and 

m)- 

Before going into details, let us motivate the use of cubature formulas in the 
present context. Let (w) denote the solution of ([T]), formally rewritten in 
Stratonovich form, if each "cii?"' is replaced by "da;*(i)", i.e. 

d 

(5) dX^icu) = (AXficu) + a{Xf{u)) - -Y^D(3,{X^{lu)) ■ P,{X^ {iu)))dt+ 

d 

i=l 

for a curve function uj = {uj-^ , . . . ,uj''') : R>o — + R'' of bounded variation. Roughly 
speaking the idea of cubature on Wiener space is to construct short-time asymp- 
totics (for some given degree of accuracy m > 2) 

N 

E{f{Xf)) = Ptf{x) = ^ifi^ti^i)) + 0{t-^), 
1=1 

with some positive, time-independent weights Ai, . . . , Xn satisfying Ai + - • ■+\n = 1 
and some d-dimensional paths wi, . . . ,ujn of bounded variation. Of course, the 
weights and paths are chosen in a specific way, which will be described later in 
more detail, and the asymptotics will only hold for some class of test functions /. 
Notice in particular that the cubature paths wi, . . . , wat depend on the interval [0, t] 
- they become rougher as t approaches 0. The aforementioned procedure replaces 
the SDE by N deterministic, well-defined PDEs, which have unique mild solutions. 
The iteration of the short-time asymptotics due to the Markov property then yields 
a weak, high order approximation scheme. 

Here also the main advantage of cubature methods in contrast to Taylor methods 
gets visible. The time-discretization in the realm of cubature methods always leads 
to reasonable expressions, namely to reasonable partial differential equations of 
type ([5]). If we wanted to apply the usual discretization methods in time like the 
Euler-Maruyama method, we might run into problems. Indeed, the naive Euler 
scheme is well-suited for the differential formulation ^ of the problem, 

t 

Xo = x and Xn = (AX„_i + a(X„_i))- + V ft(X„_i) A„B^ 

n ^-^ 

i=l 

for n > 1, but it might immediately lead to some Xn ^ T^iA). Even in the case 
of a well-defined strong solution, there is no reason why the discrete approxima- 
tion should always stay in ^{A). Hence the naive implementation of the Euler- 
Maruyama method does not work. 

Only formulation ^ seems to be suitable for using an Euler-like method. If 
one understands the semigroup S well, one can approximate X^ by expressing the 
integrals in ^ as Riemannian sums, involving evaluations of St-s, which yields 
a sort of strong Euler method (see for instance the very interesting book ^20] for 
strong convergence theorems in this direction). 

We do not discretize the integral in formulation (|3|), but (weakly) approximate 
the Brownian motion's paths by paths of bounded variation such that we obtain a 
weak approximations of X^ . 

The presence of jumps does not lead to more complicated expressions, since the 
short time asymptotics of a jump-diffusion can be easily derived from a diffusion's 
short-time asymptotics by conditioning on the jumps. The arising picture is the 
following. Discretizing the equation ^ means to allow a certain number of jumps 
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between to consecutive points in the time grid. Between two jumps we apply a 
diffusion cubature formula to express the short-time asymptotics. This yields as a 
corollary of the diffusion theory also the jump diffusion theory. 

Remark 1.1. A direct application of the cubature on Wiener space technique 
to jump diffusions driven by Levy processes of infinite activity is not possible. 
However, notice that any Levy process can be approximated by processes with 
finite activity in a weak sense, see, for instance, [3] and [19] . and therefore the 
solutions of the corresponding Hilbert-space valued SDEs converge weakly. In this 
sense, cubature methods can be use for jump diffusions driven by Levy processes 
of infinite activity, too. 

The article is organized as follows. In Section [5] we describe the analytic setting 
for a stochastic Taylor expansion to work. This is a delicate question since we 
deal with one unbounded vector field. In Section [3] and Section [4] we work out the 
cubature method from the scratch and prove the relevant convergence results in the 
diffusion case. In Section [5] we allow for jumps and prove the associated short-time 
asymptotics which is relevant to set up a weak approximation scheme. In Section 
E]we apply our method to several examples to demonstrate the results. 

2. Setting and assumptions 

Let {i},J^,P,{J^t)t>o) be a filtered probability space with a filtration {J^t)t>o 
satisfying the usual conditions. Let {Bt)t>o be a c?-dimensional Brownian motion 
and (Lj)t>o, j = 1, . . . , e be e independent compound Poisson processes given by 

fe=i 

where denotes a Poisson process with jump rate /ij > and = {Zl)k>i is an 
i.i.d. sequence of random variables with distribution Vj for j = 1, . . . , e, such that 
each i^j admits all moments. All notions are with respect to the given filtration. 
We assume furthermore that all sources of randomness are mutually independent. 

Let H he a separable Hilbert space. We furthermore fix a strongly continu- 
ous semi-group S on H with generator A. Let a, Pi, ... , fid, the diffusion vector 
fields, and 5i, . . . ,6e, the jump vector fields, be C°°-bounded on H , that is, the vec- 
tor fields are infinitely often differentiable with bounded partial derivatives of all 
proper orders n > 1. We consider the mild cadlag solution {Xf)t>o of a stochastic 
differential equation 

d e 

(6) dX^ = {AX^^ + ))dt + Y. P^^^t- )dBl -f E 5, {X^- )dLl , 

(7) X^,=xe H. 

See [5] for all necessary details on existence and uniqueness for the previous equa- 
tion. Notice furthermore the decomposition theorem in [S], which states that we 
do not need any further existence and uniqueness results in this case: in particular, 
we do not need to impose further (contractivity) conditions on A as in [B] in the 
finite activity case. 

The previous conditions are slightly more than standard for existence and unique- 
ness of mild solutions, i.e. in [HI the authors need Lipschitz conditions on the vector 
fields, whereas we assume them to be C°°-bounded. In order to formulate a sto- 
chastic Taylor expansion we shall need one main assumption, which we formulate 
in the sequel. This assumption has already been successfully applied in several 
circumstances, e.g. [2], [7] or the recent [5]. 
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We apply the following notations for Hilbert spaces 'D{A^), 

V{A'') := {he H\heV{A''-^) and A''-^heV{A)} , 

1=0 

V{A°^) = Pi V{A''), 

k>0 

which we need in order to specify the main analytic condition for our considerations: 

Assumption 2.1. We assume that a, Pi, . . . , Pd, the diffusion vector fields, and 
Si,...,6e, the jump vector fields, map 'D(A'') 'D{A^) and are C°°-bounded 
thereon for each k > 0, that is, the vector fields are infinitely often differentiable 
with bounded partial derivatives of all orders n > 1 on the Hilbert space 'D{A'') for 
each fc > 0. 

Remark 2.2. These assumptions are the appropriate analogues of the assumptions 
in finite dimension that the vector fields are bounded, C°°-bounded (see [T7] and 
|18j). In order to establish a true convergence rate one needs an additional cut-off 
argument, which is outlined in Remark l4.9l This can - like in the finite dimensional 
setting - certainly be improved. 

Example 2.3. In order to show examples of vector fields which are C°°-bounded on 
'D{A'') consider the following structure. Let _ff be a separable Hilbert space and A 
the generator of a strongly continuous semigroup. We know that is a Frechet 

space and an injective limit of the Hilbert spaces for fc > 0. Following the 

analysis as developed in ffj (see also [E] and [12] were the analytic concepts have 
been originally developed), we can consider the vector field V : U C i? — > X'(A°°). 
If V is smooth in the sense explained in [7] and has the property that its derivatives 
of order n > 1 are bounded on U C H, then V is obviously a C°°-bounded vector 
field and additionally is a Banach-map-vector field in the sense of [7]- 

Such vector fields constitute a class, where the above assumptions can be readily 
checked. 

3. The case when A is bounded linear 

We shall assume in this section that there are no jumps, i.e. we consider 

d 

(8) dX^ = {AX^ + a{Xn)dt + J2 mt)dBl. 

i=l 

In order to give an introduction to cubature on Wiener space, we consider the 
problem for a bounded operator A. In this case, there are virtually no differences 
to the finite dimensional setting, except the fact that the drift vector field does 
not need to be bounded by some constant on the whole Hilbert space (due to the 
presence of one linear operator in it). Remember that in [H] and [Hj and [I7j one 
deals with globally bounded vector fields. We shall circumvent this problem by a 
small refinement of the arguments. 

Since mild and strong solutions coincide, we can always work with strong so- 
lutions, which are semi-martingales. Consequently, we can rewrite ([5]) into its 
Stratonovich form 



(9) 



d 

dxt = p^{x-^)dt + E A(xn o dBi 

i=l 
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where Pq : H ^ H denotes the Stratonovich-corrected drift, i.e. 

1 

(10) Mx) ^Ax + a(x) - ^ ^ D/3^{x) ■ 

where 



2 

i=l 



d 

DF{x) ■v= J- F{x + ev) 

denotes the Frechet derivative of a function or vector field F. This notation enables 
us to write 

d 
1=0 

where we use the convention that ^'odB^ = dt" . 

The following notions form the core of cubature on Wiener space. We only give 
a short description and refer the reader to [18] and [22] for more details. Let A 
denote the set of all multi-indices in {0, . . . , d}. We define a degree on A by setting 

deg(ii, ...,ik)^k + if{l<j <k\ij^ 0}, 

k G N, (ii, . . . , ik) € A. We have to count all the zeros twice because of the different 
scalings for t ~ and the Brownian motion. 

Recall that any vector field /3 can be interpreted as a first order differential 
operator on test functions / by 

iPf)ix) = Dfix) ■ /3ix), xeH. 

For a multi-index {ii, . . . ,ik) £ A, k E N, let 



Br-''''= / odBll---odB,^ 

JO<ti<-<tk<t 

denote the corresponding iterated Stratonovich integral. The iterated integrals form 
the building blocks of the stochastic Taylor formula, see [T] . 

Proposition 3.1 (Stochastic Taylor expansion). Let f E C°°{H;M.) and fix < 
t < 1, m eN, m > 1, X E H . Then we have 

fixn^ (A,---A,/)(x)i?f--^^'=)-hi?™(t,/,x), 



{ii,...,ik)eA 
dog(ii,...,ifc)<m 



with 



sup yjE{Rrn{t, /, x^) < Ct 2 max sup\P^,■■■P^J{y)\. 

x£H m<dcg(ii,...,jfc)<m+2 y^H 

Remark 3.2. Notice that under the assumptions of Proposition 13 . II the bound for 
the remainder term can be infinite. However, if / is bounded, C°°-bounded and 
has bounded support, then we can guarantee that 

sup • • • A^/(y)| < oo. 
In the unbounded case this question is more subtle, see Example 14.81 and Remark 



Proposition 13.11 shows that iterated Stratonovich integrals play the role of poly- 
nomials in the stochastic Taylor expansion. Consequently, it is natural to use them 
in order to define cubature formulas. Let Cf,t,([0, i]; M'*) denote the space of paths 
of bounded variation taking values in R''. As for the Brownian motion, we append 
a component uj^{t) = t for any lo E Cto([0,t];R'^). Furthermore, we establish the 
following convention: whenever is the solution to some stochastic differential 
equation of type U]) driven by B, whether on a finite or infinite dimensional space, 
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and uj g Cbv ( [0, t] ; K*^) , we denote by (uj) the solution of the deterministic differ- 
ential equation given by formally replacing all occurrences of "odB" with ^^duj'^{s)" 
(with the same initial values), see equation ([5]) as compared to ([T]). Note that it 
is necessary that the SDE for X is formally formulated in the Stratonovich sense 
(recall that the Stratonovich formulation does not necessarily make sense for ([1])). 

Definition 3.3. Fix t > and m > 1. Positive weights Xi, . . . ,Xn summing up to 
1 and paths oji,. . . ,ujn € Cbv{[0,t];M.'^) form a cubature formula on Wiener space 
of degree m if for all multi-indices {ii, . . . ,ik) (z A with deg(ii, . . . , ik) < m, k € N, 
we have that 

N 

1=1 

where we used the convention in line with the previous one, namely 
bI''-'"'\uj) = f dLu''{ti)---duj"'{tk). 

•JO<ti<-<tk<t 

Lyons and Victoir [18^ show the existence of cubature formulas on Wiener space 
for any d and size N < #{J G A\ deg(/) < m} by applying Chakalov's theorem on 
cubature formulas and Chow's theorem for nilpotent Lie groups. Moreover, due to 
the scaling properties of Brownian motion (and its iterated Stratonovich integrals), 
i.e. 

it is sufficient to construct cubature paths for t — I. 

Assumption 3.4. Once and for all, we fix one cubature formula wi, . . . ,iON with 
weights Ai, . . . , Aat of degree m > 1 on the interval [0, 1]. By abuse of notation, 
for any t > 0, we will denote (jJi{s) — LOi{s/^/t), s G [0,t], / = 1,...,A^, the 
corresponding cubature formula for [0, t]. 

Example 3.5. For d—\ Brownian motions, a cubature formula on Wiener space 
of degree m = 3 is given by iV = 2 paths 

-i(.) = -2(.) = ^ 

for fixed time horizon t. The corresponding weights are given by Ai = A2 = ^. 

Combining the stochastic Taylor expansion, the deterministic Taylor expansion 
for solutions of ODEs on H and a cubature formula on Wiener space yields a one- 
step scheme for weak approximation of the SDE ([T]) for bounded operators A in 
precisely the same way as in 1181. Indeed, we get 



(11) sup 



£;(/(xn)-EL^'/(^*'(^')) 



<Ct 2 max sup |/3ii • • •/3i^/(?/)|, 

(ll,...,Jfc)e,4 y^H 

m<dcg(ii ,...,Zfc)<m-t-2 

for < i < 1, / e C°^{H). 

For the multi-step method, divide the interval [0, T] into p subintervals accord- 
ing to the partition Q = < ti < ■ ■ ■ < tp = T . For a multi-index (Ji, . . . ,lp) € 
{1, . . . , N}P consider the path loij^^,,,_i^ defined by concatenating the paths uji-^ , . . . , toip 
i.e. uji^,...,ipit) = uji^{t) for t e [0, ii[ and 

for r such that t G [tj.-i,ir[, where mi^ is scaled to be a cubature path on the 
interval [0,tr — tr~i\- 
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Proposition 3.6. Fix T > 0, m e N, a cuhature formula of degree m as in 
Definition \S.3\ and a partition of [0, T] as above. For every f G C°°{H) with 

sup sup • • •/3ijPt/)(?;)| < oo 

0<t<T yeH 

for all (ii, . . . ^ A with m < deg(ii, . . . ,ik) <m + 2. k (zN, there is a constant 
D independent of the partition such that 



E{f{X}))~ Az,...A,^/(X^K,...,J) 

{h,...,lp)&{l,...,N}p 



< DT max {tr - t^-i)^™"^)/^ 



sup 

xeH 



Note that the assumption on / in Proposition 13.61 is always satisfied if / is 
bounded, C°°-bounded and has bounded support and if aU vector fields a, Pi, . . . , (3d 
have bounded support (compare to Remark [3.2|) . In this case Ptf has bounded 
support, too, and all derivatives are bounded on the whole of H (see the discussion 
after Theorem 14.41 in the next section). 



Remark 3.7. Proposition 13.61 also yields determinstic a priori estimates for the 
weak rate of convergence, which hold true if we are able to evaluate the respective 
tree deterministically (see also Section E]). In principle there are methods to do so, 
see |21] . however, they might be more cumbersome to implement than a Monte- 
Carlo evaluation of the tree. 



Remark 3.8. There is one general case where Proposition 13.61 can be applied: 
namely for stochastic differential equations of type ((T|) . If we replace the unbounded 
operator A with a bounded operator A, which is close to A for a large enough set 
of values x € ^{A), then we can apply the previous result for bounded linear 
operators. One candidate for this procedure is the Yosida approximation of A. 

4. The case when A is unbounded linear 

We shall also assume in this section that there are no jumps and refer to equation 
((!]). We recall the analytic background: let us denote by 'D{A'^) the Hilbert space 
given by 'D{A'') equipped with the graph norm 



k 

|2 V^ll«,„l|2 

i=0 



h G T>{A^) and A: > 1. Furthermore, we introduce the space 

oo 

V{A°°) ^ Pi V{A''). 

'D{A°°) is topologized as projective limit of the Hilbert spaces V{A''), fc > 0, i.e. the 
topology on V{A°°) is the initial topology of the maps V{A°°) -> V{A''), fc e N. 
Note that 'D{A°°) is no longer a Hilbert space, but only a Frechet space, i.e. a 
locally convex vector space which is completely metrizable by a translation invariant 
metric. We assume our main Assumption 1 2. 1[ i.e. the vector fields restricted to the 
Sobolev spaces are C°°-bounded. 

The following proposition collects a few easy, but interesting corollaries from the 
existence and uniqueness theorem for equation ([T]) applied to the situation specified 
in Assumption 12. II 

Proposition 4.1. Fix fc e N. For any x G V{A'') there is a unique continuous 
mild solution of interpreted as an SDE in the Hilbert space D^A^). If x G 
'D{A^^^) C 'D{A^), then the mild solution in 'D{A^) coincides with the mild solution 
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in Consequently, for x £ 'D{A^'^^), the mild solution of in ^{A^) is 

a strong solution and, in particular, a semi-martingale. 

If we start in x G 'D{A°°), then we get a continuous process Xf E 'D{A°°), 
such that is the (mild and strong) solution of ([5]) in any Hilbert space I?(A'^), 
k € N. Furthermore, Proposition 14.11 allows us to avoid any problems due to 
the topological structure of 'D{A°°) by reverting to an appropriate Hilbert space 
'D{A'') and interpreting the results again in 'D{A°°). The meaning of I?(yl°°) is that 
it is the largest subspace of the Hilbert space H, where we can innocently do the 
necessary analysis on differential operators. Notice that there are subtle phenomena 
of explosion, which can occur in this setting: for instance the law of a strong solution 
process X solving equation ([T]) might be bounded in H but unbounded in TJ(A), 
where it is a mild solution. Due to such phenomena. Example 14.81 and Remark 14.91 
after Theorem 14.41 are in fact quite subtle. 

As in Section [HI we introduce the vector field /3o defined by 

1 

(12) Po{x) =Ax + a{x) - - ^ Dp,{x) ■ /3,{x). 

i=l 

Po is defined for x € T^iA). As a vector field taking values in 'D{A'^), it is only 
well-defined on V{A''+^). Consequently, for x £ V{A''+^), we may reformulate the 
SDE dl]) - understood as equation in 'D{A'') - in Stratonovich form 

Now we formulate the stochastic Taylor expansion in some I?(A'"^™^) with a 
degree of regularity r(m) > depending on m > 1. For the estimation of the error 
term, we will use the extended support esupp(Xj^; loi, . . . , lu^) defined by 

(13) esupp(Xf ; wi, ...,u;n)^ supp(Xn U {X^{lui), . . . (c^jv)}, 

where t > 0, x £ H , and toi, . . . ,ujn are paths of bounded variation. Here supp(Xf ) 
means the support of the law of in I?(A'''™^). Despite Assumption 13.41 let us, 
for one moment, enter the dependence of the cubature formula on the interval [0, t] 
explicitly into the notation, in the sense that w^*'' , • . . , w^'' are the paths of bounded 
variation scaled in such a way that they, together with the weights, form a cubature 
formula on [0,t]. Then we denote 

6t(x)= y esupp(X,-;c.i*\...,4*)). 

0<s<t<T 

Remark 4.2. If a general support theorem holds in infinite dimensions, we can 
replace the extended support by the ordinary support of , since the solution of 
the corresponding ODE driven by paths of bounded variation lie in the support 
of the solution of the SDE according to the support theorem. However, up to our 
knowledge, no general support theorem has been established for our setting so far. 

Theorem 4.3. Let m > 1 be fixed, then there is r{m) > such that for any 
f £ C°°{H;R), X £ X>(A''(™)), and < t < 1 we have 

/(Xf)= J2 {l3,,---l3,J){x)Bf'-'"'HR^{t,f,x), x£V{A^^'"% 

deg(ii,...,ifc)<m 

with 

VE{RJtl^) < Ct"^ max sup |i?(A, • • • AJ(X,^))| 

m<dcg(ji,...,ifc)<m+2 0<s<t 

< Ci"^ max sup \I3^, ■ ■ ■ I3^J {y)\ . 

m<dcg(ii,...,ife)<m+2 j;gSt(2^) 

We can choose r{m) = + 1, where [yJ is the largest integer smaller than 
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Proof. The proof is the same as in the finite dimensional situation, but one has to 
switch between different spaces on the way. 

Fix m and / as above and x £ I?(aL"2"J+^). We interpret the equation in 

2"J+^). By the above remarks, we can express the SDE in its Stratonovich 
form By Ito's formula, 

(14) f{Xn = fix) + [\l3of){X^)ds + E fm){Xf) ° dBl. 

Jo Jo 

The idea is to express {/3if){X^) again by Ito's formula and insert it in equa- 
tion (fH|) . This is completely unproblematic for i G {1, ■ • ■ ,d}. For i — 0, recall 
that 

1 

(/3o/)(x) = Dfix) ■ Ax + Dfix) ■ (a(x) - - ^ Df3^ix) ■ A (a:)) . 

i=l 

By re-expressing ()14p in Ito formulation, applying Ito's formula, and re-expressing 
it back to Stratonovich formulation, we see that 

{Pof)ixf) = (/3o/)(x) + (\[3if){x:)du + E fmof){x:) o dsi, 

Jo Jo 

where 

(15) {f3lf ){x) = D^f{x){Ax, Ax) + Df(x) ■ {A^x + Aa{x) + ■■■) + ■■■, 

provided that all the new vector-fields are well-defined and the processes {f3i(3of){Xi^) 
are still semi-martingales. Both conditions are satisfied ii x £ 'D{A'^) - notice that 
the maps 

p(-^fc+i) _^ V{A''), X ^ Ax a,re C°°, k e N. By induction, we finally get 

AJ)(x)S^-'^'+i?„(t, f,x) 



(n,---,ik)eA,t„e{o,...,d} Jo<to<-<tk<t 

dog(n,...,ife)<m<dGg(io,ii,...,ifc) 

Note that Rm is well-defined for x G I?(y4L'2"J ''"^) because integration of non-scmi- 
martingales with respect to dt is possible, which corresponds to the index ig = 0. 

As in the finite dimensional case, we re-express Rm in terms of Ito integrals and 
use the (one-dimensional) Ito isometry several times, until we arrive at the desired 
estimate. □ 

We recall the notation Ptf{x) = E{f{Xf)) for bounded measurable functions 
f : H ^ M.. Analogously to Proposition 13.61 we immediately get the following 
theorem. 

Theorem 4.4. Fix T > Q, m > I, r{m), x G X>(A''('")) as in Theorem\Jl^ a 

cubature formula on Wiener space of degree m as in Definition \3.3\ and a partition 
= to < ti < ■ ■ ■ < tp = T . Under Assumvtion \2.1\ for any f G C°°{H; R) with 

sup sup • • •/3,^Pt/(y)| < GO 

0<t<T yeerix) 



f{xf) ^ J2 (^^1 



{ii,...,ik)eA 
eg(ii,...,ifc)<m 



with 

Rm{t,X,f) 



E 
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for all (ii, . . . ,ik) € A with m < deg(ii, . . . ,ik) <m + 2. A; G N, there is a constant 
D independent of the partition such that 

(h,...,lp)(i{l,...,N}p 

<DT max (t,. - t^_i)(™-i)/2, 

r—l,...,p 

where X^{u)) is, again, understood as the mild solution to an ODE in for 
any path uj of bounded variation. 

Proof The proof follows Kusuoka [16] , , see also [M] . For / : — > M and x e H 
let 

N 

Q(t)f{x) = Y.\if{X^{uJi)), 

1=1 

where uji, . . . ,lui are scaled to form a cubature formula on [0,t]. Denote At^ = 
tr — tr-i, f — ■ ■ ■ iPi the increments of the time partition given in the statement 
of the theorem. By iterating the operators Q(Atr) (and the semigroup property of 
ODEs), we immediately obtain 

(16) E K---\f{x^i^h.....,0)^QiAt,)°---°QiAt,)f{x). 

(/i,...,ip)e{l,...,Ar}p 

By ordinary Taylor expansion, keeping in mind the degree function deg, we note 
that 

N 

Q(t)f{x)^Y.^l E {P^,■■■P^J){x)Bf'-''''\u;l)+Rm{t,X,f), 

deg(ii,...,ifc)<m 

where x G 2?(A''(™)) and 

N 

R^(t,x,f) = Y.^i E 

1=1 (ii,...,ifc)e^,ioe{0,...,d} 

deg(zi,...,ifc)<m<deg(io,---,jfc) 

/ (A„ • ■■P,J){Xl{ui))dw\''{t^) ■ ■■d^l^itk). 

JO<to<-<tk<t 

In the sequel C denotes a constant independent of the partition and x, which may 
change from line to line. We can estimate the approximation error by 

<C sup max |A„ • • • ft J(Xf (^O)I ■ 

0<s<t, 1 = 1, ...,N m<dog(jo,....ifc)<m+2 

Combining this result with Theorem 14.31 we may conclude that 

(17) \Ptfix)~Q^t)fix)\<C sup \P,„---P,J{y)\t'^. 

By telescopic sums, 

Prfix) - (3(Atp) o • • • o Q(Ati)f{x) = 
p 

E Qi'^tp) ° ■ ■ ■ ° Q{At,.+ l)iPtrf{x) - Q(At,^)Pu^J{x)). 

r=l 

For the estimation of the rear term 

PtJ{x) - QiAt^)Pu.J{x) = (Fa*. - Q^Au))Pu^^f{x), 



Rni{t,X,f) 
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we may use PTl) with f{x) being replaced by Pt^_^f{x), giving us 

p 

\PTf{x)~Q^At,) 0---0 Q^At,)fix)\ < J2 ^''^ \Puf{x) - QiAt^)Pu-^f{x)\ 

r=l 

P 

<CY^ sup \(3.„---p,,Pt,._J{y)\iAtr)'^ 

<c sup |Ao---A.J^t/(j/)lE(Ai.)'^, 

yeeT{x),o<t<T 

from which we may easily conclude the theorem. □ 

Remark 4.5. Gyongy and Shmatkov [TT] show a strong Wong-Zakai-type approx- 
imation result, where they also need to impose smoothness assumptions on the 
initial value x. Otherwise, the assumptions in [11 are different from ours. They 
allow linear, densely defined vector fields and general adapted coefficients, on the 
other hand the generator A needs to be elliptic. 

Remark 4.6. Under the previous assumptions we can also prove a Donsker-type 
result on the weak convergence of the "cubature tree" to the diffusion. This result 
will be presented elsewhere. 

Remark 4.7. If / is smooth then we can show by (first and higher) variation 
processes, as introduced for instance in [8], that x i-^ Ptf{x) is smooth on 'D{A^). 

Fix > 0. Let Jo^t{x) ■ h denote the first variation process of in direction 
h G V{A''), i.e. 

£=0 

Jo-,t{x) ■ h is the mild solution to an SDE of the type ([8]). Consequently, it is 
bounded in L^{il,, !F, P;'D{A'')) and we may conclude that 

|- Ptfix + €h) = E(Df{Xn ■ Jo^tix) ■ h) 

exists and is bounded by boundedness of Df and intcgrability of the first variation. 
Similarly, we get existence and continuity of higher order derivatives on 'D{A'^). 



Jo^t{x) ■ h 



Example 4.8. We shall provide examples, where the assumptions of Theorem[ 
are satisfied, i.e. where we obtain high-order convergence of the respective cubature 
methods. The conditions seem at first sight restrictive (see the following Remark 
14.91 for a concrete example under Assumptions [2T|) . however, the conditions are 
parallel to those obtained in and [Hj, where the functions and vector fields 
have to be bounded and C°°-bounded. Here we have an additional complication of 
one certainly unbounded, but not even continuous drift vector field, which leads to 
the following set of assumptions. 

The vector fields a, (3i, . . . , Pd have the following property (compare also to tame 
maps in '12]): there exist smooth maps 

5,^l,...,^rf:I?(^''(™))^I?(A2r(m)) 

such that for A: > the restrictions to the respective subspaces I?(A'^+''(™)) C 
5, A, • ■ • : I?(A'=+''(™)) ^ p(^'=+2'-(™))^ 

and such that 

(18) a^ao (i?(A, A)'-^")), A = A ° (i?(A, A)'^^")), 
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for z = 1, . . . , d. Here i?(A, A) denotes the resolvent map for A S p{A). We assume 
that a, (3i, . . . , (3d have bounded support on The function / is of the 

same type f — g o (i?(A, A)'^'-™-') for a bounded, C°°-bounded function g : iJ ^ R. 

Under these assumptions we can readily check that the law of the mild solution 
Xf starting at the initial value x £ 'D{A^^™'>) has bounded support in H: outside a 
ball of radius R > in H the solution process becomes deterministic, = Stx on 
some interval, hence by the uniform boundedness theorem there is a large R' such 
that the image of the ball with radius R > under the maps St lies in a ball with 
radius iJ' > on [0,T]. 

For smooth functions / of the stated type we then have 

sup sup • • • A,Pt/(2/)| < oo. 

0<t<T yeH,\\y\\<R' 

Since we only take the suprcmum over bounded sets, namely the extended supports 
of , this implies the assumption of the Theorem 14.41 

Remark 4.9. The previous assumptions (fT5|) on the vector fields are not too re- 
strictive since we can always obtain them by a linear isomorphism and (smoothly) 
cutting off outside a ball in 'D{A^'^"^^). Both operations are numerically innocent. 
Under Assumption 12.11 we can apply the following isomorphism to the solution of 
our SDE (dl: 

This isomorphism transforms the solution on 'D{A^^"^'') to an iJ- valued process 

y,^ -i?(A,A)-'^(")xf(^"^^'^'"'^ 

with X = R{X, AY^"^^y. The transformed process satisfies an SDE, where the 
transformed vector fields (if well defined) factor over i?(A, AY^"^^ such as previously 
assumed in the assumption (jlSp . namely 

(19) dyy = {Ayy + {{r{\, a)""^"^ o a o r{\ Ay'-"'^){Yyt))dt+ 

d 

+ Y.{R{\ A)-'^("' o 13, o i?(A, A)'-("))(r/)dB^ 

i=l 

The assumptions (IT5|) mean that we must (smoothly) cut off the vector fields 
a, /?!,..., /3d outside sets of large norm |Mlu(^r(„j')'), which is an event - under 
Assumption 12.11 - of small probability (recall that the vector fields are Lipschitz on 
2?(^''("i)) and therefore second moments with respect to the norm |Mli5(^r(„j)) ex- 
ist). Notice that the cut-off vector fields do not have an extension to H since contin- 
uous functions with bounded support on I?(v4''(™)) do generically not have a contin- 
uous extension on H . For {Y^ )^y^^ we can take initial values y Cz H, however, those 

initial values correspond to quite regular initial values x — R{X, A)'^''™^y e V^A^^"^^) 
for the original process (-^f )t>o- 

From the point of view of the process Y we have hence proved that 

/(r#) = go (i?(A, A)'^('")) o (i?(A, A)-^("))(X|;) = g(X^) 

is weakly approximated by evaluating / on the cubature tree for Y. This is equiva- 
lent to evaluating g on the cubature tree for X in order to approximate the expected 
value E{giX^)). 



14 



CHRISTIAN BAYER AND JOSEF TEICHMANN 



5. The cubature method in the presence of jumps 

The extension of cubature formulas to jump diffusions seems to be new even in 
the finite dimensional case. We shall heavily use the fact that only finitely many 
jumps occur in compact time intervals almost surely. 

We shall first prove an asymptotic result on jump-diffusions. 

Proposition 5.1. Consider equation Let f : H W be a bounded measurable 
function, then we obtain 



(20) E{f{Xn)^ J2 ^ll-^e-*^^"- 

— ^ nil ■ • • Tie! 
ni ,. ..,ne >0 

xE{f{Xn\Nl=n, forj^l,...,e) 

for t > 0. 

Proof. We condition on the jump times and read off the results by inserting the 
probabilities for a Poisson process with intensity fij to reach level Uj at time t. □ 

This result gives us the time-asymptotics with respect to the jump-structure. It 
can now be combined with the original cubature result for the diffusion between 
the jumps, in order to obtain a result for jump-diffusions. We denote by the 
jump-time of the Poisson process for the n-th jump. We know that for each 
Poisson process the vector (rj* , . . . , t^^. — t;^^._i, t — r-j^. ) is uniformly distributed if 

conditioned on the event that = rij > 1. The uniform distribution is on the 
rij-simplex iA"j C This allows us to apply an original cubature formula 

between two jumps of order m — 2ni — ... — 2ne, since we gain, for each jump, one 
order of time-asymptotics from the jump structure. 

Assume now that the jump distributions Vj are concentrated at one point Zj ^ 0, 
i.e. AL-'j — Zj for j — l,...,e and k > 1. If we want to consider a general 

jump-structure this amounts to an additional integration with respect to the jump 
distribution Vj. 

Now we define an short-time approximation for the conditional expectations 

E{fiXn\Ni=n, forj = I,...,e) 

of order m — 2ni — ... — 2ne with ni + . . . + < "^2^ . Expressed in words, 
we are going to do the following: starting from the initial value x € 'D{A^^"^^) 
we solve the stochastic differential equation ([2]) along the cubature paths ui with 
probability Aj > 0, I — 1,...,N until the first jump appears. We collect the 
end-points of the trajectories, add the jump size at these points and start a new 
cubature method from the resulting points on. Notice that we can take a cubature 
method of considerably lower degree since every jump increases the local order of 
time-asymptotics by I. The jump times are chosen independent and uniformly 
distributed on simplices of certain dimension Uj such that ni + . . . + rte = n < 
We denote the cubature trajectory between jump t^_-^ and r| for I < q < Uj with 
oJij^q. If nj = no trajectories are associated. Hence we obtain the following 
theorem: 

Theorem 5.2. Fix m > I. Consider the stochastic differential equation 0) under 
the condition — Uj for j — 1, . . . , m with ni + . . . + n,, ^ n along concatenated 
trajectories of type toij^q- Choose a cubature method of degree 

m' — m ~ 2n > 1. 

Concatenation is only performed with increasing q-index and a typical concatenated 
trajectory is denoted by uji-^^...^i^. Here we have in mind that the intervals, where 
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the chosen path is (^i.j,q, come from a jump of and have length — t^^-^^. Then 
there is r{m') > such that 

N 

(21) \E{f{Xn\Ni^n,)- ^■■■^i.E{f{X!{u^i,,..jJ)\Ni =n,)\ 

ii,...,i„=i 

I^Ct"^ max sup • • • A,£;(/(X^^,,(^i,^„...,i„)) | TVf = n,)|, 

{ii,...,ik)€A yGsupp(XJ), 
dcg(ii,...,ifc)<m+2 0<s<t 

where (w) means the solution of the stochastic differential equation Q in Strato- 
novich form 

d e 

(22) dX^iu;) = f3o{X^- {u^))dt + J] A(Xf_ {u;)) o dBl + ^ S,{X^^ )dLl 

i=l 3=1 

along the trajectory lu. 

Proof. By our main Assumption 12.11 we know that the linkage operators x i— > (x) 
are C°°-bounded on each 'D{A''), hence through concatenation the errors, which 
appear on each subinterval [tn-q,t[ are of the type y i?(X^^_j(a;;^_^jg_...^i^)) for 
some 1 < q < n. Taking the supremum yields the result. □ 

Combining the previous result with Proposition 15.11 yields under certain con- 
ditions on the vector fields (see the discussion after Theorem 14.31 in the previous 
section) by the triangle inequality that there is a constant D > such that 

N 

(23) \E(f(xn)~ y y -t^ini-...-t^,n,^ 

2(ni + ...+ne)<m ^' 

xAi, . . . Ai„i;(/(Xf (L^z,,...,iJ) I Nl = n,)\ < Dt"^ < oo, 

By iteration of the previous result we obtain in precisely the same manner as in 
Section 2] a cubaturc method of order m by applying several cubature methods of 
order m' < m between the jumps. 



Remark 5.3. The only random element in the expectation £^(/(Xf^(a;;j^...^;^)) | 
%j) is given by the jump times r^, 



rij) is given by the jump times r^, which vary on certain simplices. For the imple- 



mentation one has to simulate the uniform distributions on the simplices iA'^'. Since 
the integrals on the simplices i A*^' only have continuous integrands, we cannot hope 
for other methods than Monte-Carlo. The evaluation by a Monte-Carlo-algorithm 
can also be seen as a random choice of the concatenation grid for the constructed 
t^ii, ...,;„• Another view could be to see a deterministic grid for the diffusion which 
is saturated by points where jumps occur. Implementation will be done elsewhere. 

6. Numerical examples 

Due to the previous results a numerical scheme for equation ([T]) can be set up by 
the following steps. Notice that we have a weak order of approximation ^^^^^ only 
under the assumptions of the previous section. In order to obtain those assumptions 
one has to modify a general equation of type (jT]) by smoothing procedures. These 
modifications can be done in a controlled way, more precisely, for each modification 
we have a rate of convergence to the un-modified object. 

• Approximate the vector fields a, (3i, . . . , Pd by vector fields satisfying the 
Assumptions 12.11 If the original vector fields are globally Lipschitz one can 
do this approximation with a rate of convergence for the L^-distance of the 
original solution process and its approximation. 
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• Choose a degree of accuracy m > 2, which determines the weak order of 
convergence in the sequel. Associated to m the number r{m) can be 
identified, which teUs us about the degree of regularity of the initial value 
X G H, which one needs for the assertions of Theorem 14.41 

• Identify due to the previous specifications a radius R > such that the 
■D(A''('"))-norm of is rarely beyond R. Cut-off the vector fields smoothly 
on and verify assumptions (fT8|) - maybe after smoothing - for 
the transformed process Y such as exercised in Remark 14.91 

• The resulting tree of trajectories yields a finite number of non-autonomous 
PDEs (O, which have to be evaluated. In the implementation one calculates 
with the smoothened vector fields satisfying the Assumptions 12.11 but not 
with the cut-off vector fields, since a large I?(A''(™))-norm value R is reached 
with small probability due to its very choice. 

We test the cubature method for two concrete examples: one toy example, where 
explicit solutions of the SPDE are readily available, and another, more interesting 
but still very easy example. Since cubature on Wiener space is a weak method, we 
calculate the expected value of a functional of the solution to the SPDE in both 
cases, i.e. the outputs of our computations are real numbers. 

The results presented here are calculated in MATLAB using the built-in PDE- 
solver pdepe for solving the deterministic PDEs given by inserting the cubature 
paths into the SPDE under consideration. This PDE solver depends on a space 
grid given by the user as well as on a time grid, which is not very critical because 
it is adaptively refined by the program. 

We do not use recombination techniques for cubature on Wiener space as in [3T] 
and use the simplest possible cubature formula for d = I Brownian motions: 

with weights Xi — X2 ~ define a cubature formula of degree m — 3 on [0, T]. 
Consequently, solving an SDE on a Hilbert space with cubature on Wiener space for 
the above cubature formula and p iterations means solving 2^ PDEs. This starts to 
get restrictive even for a very simple problem for, say, p = 10 - where one already 
has to solve more than one thousand PDEs. One possibility to overcome these tight 
limitations is to use "a Monte-Carlo simulation on the tree" . Recall that an p-step 
cubature method approximates 

miX^)) « E ■ ■ • A,,/(X^K,...,,,)). 

Since ^ A^^ • • • Xj^ = 1, we can interpret the right hand side as the expectation of 
a random variable f{X^{uj.)) on the tree {1, . . . , N}p. Therefore, we can approx- 
imate the right hand side by picking tuples (ji, . . . ,jp) € {1, . . • ,N}p at random 
- according to their probabilities Xj^ ■ ■ ■ Xj^ - and calculating the average of the 
corresponding outcomes f{X'^{LUj^^,,,^j^)). Of course, by following this strategy we 
have to replace the deterministic error estimates by a stochastic one, which heavily 
depends on the standard deviation of J{X^{uj.)) understood as a random variable 
on the tree. Notice, however, that this is the usual situation for weak approximation 
methods for SDEs. 

Consider the Ornstcin-Uhlenbeck process Xf defined as solution to the equation 
(24) dXf = AATf -)- <l)dBt 

on the Hilbert space H = i^(]0, A denotes the Dirichlet Laplacian on ]0, 1[, 
i.e. A is a negative definite self-adjoint operator on H with I?(A) — Hq{\Q, 1[) n 
i/^(]0, 1[) extending the classical Laplace operator defined on C^(]0, It is easy 
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to see that A is dissipative and therefore, by the Lumer-PhiUips theorem, it is the 
generator of a Co contraction semigroup {St)t>o on H. The coefficient (j) G H is 
some fixed vector. 

fn this case, the definition of a mild solution 

(25) = Stx+ f St-sHBs 

Jo 

already gives a representation of the solution provided that the heat-semigroup St 
applied to the starting vector x and to (p is available. We choose x{u) = sin(7ru), 
u g]0, f [, and may conclude that 

Stx = e-'^^^x 

because x is an eigenvector of A with eigenvalue —tt"^. Consider the linear functional 
^ : H given by 

(26) $(y) = / y{u)du, y G H. 



We want to compute 

E{<^{Xf)) = E(^j e-''\m{TTu)du + J J Si^s<l>{u)dBsdu 

= / e"''' sm{TTu)du = 0.3293 x 10"^ 
Jo 



p 


Error 


1 


-0.3601 X 


10- 


-4 


2 


-0.2192 X 


10" 


-4 


3 


-0.1226 X 


10- 


-4 


4 


-0.0652 X 


10- 


-4 


5 


-0.0334 X 


10- 


-4 


6 


-0.0172 X 


10- 


-4 


7 


-0.0084 X 


10- 


-4 


8 


-0.0031 X 


10- 


-4 


9 


-0.0002 X 


10- 


-4 


10 


-0.0013 X 


10- 


-4 



Table 1. Error for the cubature method in the OU-case (absolute error) 



In Table [1] the error, i.e. the output of the method minus the true value given 
above, is presented, p is the number of cubature steps, i.e. the number of iterations 
of the one-step cubature method. The discretization in space, i.e. of ]0, 1[, used by 
the PDE solver contains 50 uniform points, the discretization of the time-interval ~ 
additional to the one induced by the cubature method - contains 500 points. The 
stochastic perturbation factor ip is chosen to be </)(«) = sin(7ru), i.e. (j) S 'D{A°°) 
even. We see a very fast decrease of the error in this simple situation. On the other 
hand, the variance of the random variable on the tree considered before is clearly 
too high for the Monte-Carlo simulation on the tree to work. Indeed, $(Xf ) has 
true standard deviation of 

(27) sd(<I>(Xf )) ^ yjj(l-e-2-^) - 0.1433. 

Assuming that the central limit theorem applies, confidence intervals around the 
solution given by a Monte-Carlo method are proportional to the standard deviation 
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Figure 1. Relative errors for the OU-process (equatfon ([M])) and 
the process with nonUnear volatiHty (equation (I28p). In the latter 
case, both initial values x{u) = sin(7rM) - referred to as "regular 
case" - and x given by ((5T|) - referred to as "irregular case" - are 
used. 



divided by the square root of the number of trajectories. Consequently, we would 
roughly need to calculate 10^^ paths on the tree in order to achieve a similar level 
of exactness as in Table [T] Indeed, note that the standard deviation of the solution 
is of order 10~^, while the error in the last row of Tabled is of order lO"'^. The 
equation 




then gives m ~ 10^^. Note that this heuristics is also confirmed by our experiments, 
where Monte-Carlo simulation on the tree clearly fails. The data are also shown in 
Figure [H 

Remark 6.1. The failure of Monte-Carlo simulation on the tree also applies to 
any other (naive) Monte-Carlo approach to problem , including the usual finite 
element or finite difference approaches. 

As a more realistic example we consider the heat equation with a stochastic 
perturbation involving a Ncmicky operator. More precisely, consider 

(28) dX^ = AX^dt + {smoX^)dBt, 

with x{u) = sin(7ru). Even though we do not know the law of the solution Xf 
of we are still able to calculate E{^{Xf)) explicitly because $ is a linear 

functional. Indeed, Xf is given by 

(29) Xf = Six+[ Si-s{sinoX^)dBs 

Jo 



CUBATURE ON WIENER SPACE IN INFINITE DIMENSION 



19 



and, consequently, 

= <t>{Six) + [ <t>{Si^,sinoX^)dB,. 
Jo 

The expectation of the (one-dimcnsional) Ito-integral is and we get the same 
result as before, i.e. 

E{<P{X^)) = <i>{Six) = 0.3293 x lO"" 

for x{u) = sin(7rM). Nevertheless, we believe that this example is already quite diffi- 
cult, especially since the cubature method actually has to work with the Stratonovich 
formulation 

(30) dX'^{uj) = (AX^iuj) - ^{cosoX^{uj)){smoX^{Lu)))dt + {sin o X^ {Lj))dLu{t). 

In particular, the equation (in Stratonovich form) has a non-linear drift and a 
non-linear volatility. 

Note that we expect the standard deviation of the solution of the above equation 
to be smaller than before, because (sinoXf)^ decreases as decreases in t. The 



I 


Error 


1 


-0.2907 X 10- 


-4 


2 


-0.2163 X 10" 


-4 


3 


-0.1467 X 10" 


-4 


4 


-0.0961 X 10" 


-4 


5 


-0.0622 X 10" 


-4 


6 


-0.0385 X 10" 


-4 


7 


-0.0228 X 10" 


-4 


8 


-0.0142 X 10" 


-4 


9 


-0.0086 X 10" 


-4 


10 


-0.0040 X 10" 


-4 



Table 2. Results of the cubature method for pS)) (absolute error) 

space discretization used by the PDE-solver has size 50, which already seems to 
be sufficient, because using a finer discretization (100 grid points) does not change 
the results significantly. Table [3] shows the results using Monte-Carlo simulation 



I 


m 


Error 


Stat. Error 


5 


32 


0.0567 X 10" 


-4 


0.1498 X 10" 


4 


10 


1000 


-0.0325 X 10" 


-4 


0.0179 X 10- 


4 


15 


1500 


-0.0184 X 10- 


-4 


0.0172 X 10- 


4 


20 


2000 


0.0128 X 10- 


4 


0.0170 X 10- 


4 


25 


2500 


0.0179 X 10- 


4 


0.0145 X 10- 


4 


30 


3000 


0.0596 X 10- 


4 


0.0167 X 10- 


4 



Table 3. Results of the cubature method with Monte-Carlo sim- 
ulation on the tree for ([28|) (absolute error) 

on the tree, m denotes the number of trajectories followed, while the "Statistical 
Error" in the table is an indicator for the error of the Monte-Carlo simulation. 
More precisely, the values in the last column are the empirical standard deviations 
of the result divided by the square root of the number of trajectories. Comparable 
to the Ornstein-Uhlenbeck process, the convergence of the pure cubature method 
is very fast, see Table[5J The (empirical) variance is, however, quite large such that 
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the Monte-Carlo aided method does not work at all. Note that the statistical error 
in Table [3] is of the order of the total computational error, which can be almost 
completely attributed to the Monte Carlo simulation. 

To test the method further we also try more irregular data. Let 



(31) 



x{u) = 



1 



1 



The exact value of the quantity of interest E{^{Xf)) = $(S'ix) is calculated 
by solving the corresponding heat equation numerically. This gives the value 
E{^{Xf)) = 0.3002 X lO-'^. The initial vector x given in ^ is in i^(]0, 1[) 
but its derivative is no longer square-integrable. Consequently, x ^ 'D{A) and the 
theory does not provide an order of approximation. Nevertheless, probably due to 
the smoothing-properties of the Laplace operator numerical results show the same 
behavior as before, see Figure [T] 

If we replace the heat equation by an evolution equation of the form 



(32) 



dXf = —Xfdt 
du 



sino XfdBt 



then we still see the same kind of behavior if we fix the space-discretization for the 
PDE-solver. This time, the PDEs require a much finer space resolution in order to 
give reliable numbers. 
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